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I. INTRODUCTION 

While it has been shown that wave maps on a 1 + 

1 dimensional Minkowski spacetime base evolved from 
smooth initial data exist for all time , and that those 
on an m + 1 (m > 3) Minkowski spacetime base can 
blow up in finite time |^ , global existence for the 2 + 1 
case remains as yet miresolved. Scaling considerations 
identify 2 + 1 as the critical dimension for wave maps, and 
so there is considerable interest in determining if indeed 

2 + 1 wave maps developed from smooth initial data can 
become singular in finite time or not. Here, we describe 
numerical work which strongly supports the contention 
that, at least for some sets of smooth initial data, they 
can. 

There are special classes of 2 + 1 wave maps for which 
global existence has been shown to hold: a) spherically 
equivariant wave maps with convex [Q, or slightly more 
general targets b) spherically symmetric wave maps 
with compact targets (plus a further technical condition 
on the target) Q, c) general wave maps (general target) 
with sufficiently small energy. 

Not included in any of these three classes are spheri- 
cally equivariant wave maps from 2 + 1 Minkowski space- 
time into the round two-sphere with initial data of ar- 
bitrary energy. Shatah and Struwe Q have conjectured 
that singular behavior should be found in this class. Our 
numerical results reported here strongly support the va- 
lidity of this conjecture. 

We examine one-parameter families of data, with small 
values of the parameter corresponding to small energy 
data and therefore global existence, and with large values 
of the parameter corresponding to data possibly leading 
to singularity formation. One might hope to find espe- 
cially interesting wave map development for data at or 
near the transition between small and large values. While 
this sort of "critical" behavior has been seen and studied 
in 3 + 1 wave maps ||^,^ , we have not found nearly as 
clear an indication of universal critical behavior for the 
present 2 + 1 case. This criticality issue needs further 
study, and is not treated in this paper. Here, our focus 



is on numerical evidence for singular wave map evolution 
from regular initial data. 

We note that our studies of singularity formation in 2+ 
1 wave maps have been carried out independently of the 
work of Bizoh, Chmaj, and Tabor |l0|| using numerical 
algorithms which differ from theirs. However their results 
and ours agree substantially. 

II. THE EQUATIONS 

Generally a wave map is defined to be a map cf)"^ from a 
spacetime (the "base" ) into a Riemannian geometry (the 
"target" ) , with (j)"^ a critical point for the action 

S[^] - / v'^^'gAB (0) {d^cl^^dA'') (1) 

where qab is the Riemannian metric on the target mani- 
fold iV" , and r]^'^ is the (inverse) Lorentz-signature metric 
on the spacetime M™"*"^. The Euler-Lagrange equations 
for this action take the form 

5^a^(/.^ + r^caM'/'''9'''/''' = o (2) 

where represents the Christoffel symbols corre- 

sponding to the target metric qab- This is a semilinear 
hyperbolic PDE system for 0^. We note that for certain 
targets, wave maps are known to physicists as "nonlinear 
sigma models." 

As noted above, the case of primary interest here is 
2+1 Minkowski spacetime for the base and the round 
two sphere for the target. In this case, the wave map 
PDE system (H) may be rewritten in the following form 

□ (/.•^ + (9^/9^0^) Sbcc^" = 0. (3) 

where the indices a, 6, c take the values {1, 2, 3} (indexing 
the ambient Euclidean 3-space for the target two sphere), 
and 5bc is the metric for this ambient space. If we now 
impose the condition that the maps be spherically 
equivariant with angular wrapping number k, and write 
(ff'{r,9,t) in the "hedgehog" form 
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sinx(r, t) sin kO 

sin xir,t) cos k9 | (4) 
cos x{r,t) 

where r is the radial distance from the origin and 6 
the azimuthal angle, then the wave map PDE system 
reduces to the single equation 



1 fc2sin2x 
X=-{rx) 



(5) 



where a prime and a dot denote partial derivatives with 
respect to r and t respectively. Thus the study of the 
Cauchy problem for 2+1 spherically equi variant (fc- 
wrapped) wave maps into the round two sphere focuses 
on finding solutions xl'^iO ^o Eq. with regular ini- 
tial data x(?', 0), x(^jO)- Note that regularity at r = 
requires that we set x(0, t) = for all t. 

While it may be interesting to examine if there is any 
variation of the behavior of solutions for wrapping num- 
bers k greater than one, we restrict our attention here to 
the single angular wrapping case k = 1. 

As for any field theory on Minkowski space, there is a 
divergence- free stress-energy tensor T^^ associated with 
wave maps. From T^,y, we obtain the energy density func- 
tion for spherically equivariant wave maps 



P{r,t) 
whose integral 



sin^X 
2r2 



E{t) 



p{r, t) r dr 



(6) 



(7) 



is conserved (ie, E{t) = E{0) for all t). The energy is a 
useful monitor of numerical accuracy, as discussed below. 



III. NUMERICAL STUDIES OF SINGULARITY 
FORMATION 

Our numerical experiments consist of specifying 
parametrized families of initial data {xa(?', 0), xaI?", 0)} 
and numerically evolving a number of sets of such data 
in each family. A typical family-one of the simplest-is 
the approximately ingoing Gaussian pulse 



X(r,0) 



x'(^o). 



(8) 



This family has three parameters A, i?o, and 6, with the 
most important one for our discussion being the scale 
parameter A. Note that the ingoing character of these 
solutions, which results from the choice of x{r,0), mini- 
mizes outer boundary effects. Note also that while, ana- 
lytically, x(0, 0) is not zero, for the choices of Rq and 5 
which we make, we can force x(0, 0) to be zero and retain 
smoothness to within numerical accuracy. 



We evolve using a second order finite difference approx- 
imation to Eq. (^) . We use an iterative Crank-Nicholson 
scheme implemented with RNPL and also make 

use of the adaptive mesh framework developed by Chop- 
tuik We have verified that the code generates so- 

lutions which converge quadratically in the grid spacing 
and conserve energy. In arguing that we are indeed gen- 
erating singularities, we will discuss the convergence and 
energy conservation tests in more detail below. 
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FIG. 1. Snapshots of the energy densities (times r) for 
a single family of initial data with varying amplitude. At 
t — the initial energy densities corresponding to ingoing 
Gaussian initial data {Ro — 8, 5 — 2.3) are shown. Super- 
critical {A = 1.4) is shown long-dashed (red) which exists 
only until t = 6. Sub-critical {A = 1.0) is shown solid (ma- 
genta). Two near-critical evolutions are also shown: slightly 
subcritical {A — 1.19) shown dotted (blue) and slightly su- 
percritical {A = 1.195) shown short-dashed (green). The two 
near-critical evolutions coincide at the scale of this graph until 
t = 12 after which we cannot compute the apparently singu- 
lar supercritical solution. The energy densities reached by the 
supercritical solutions extend significantly off the scale of this 
graph. 

For a general set of ingoing Gaussian pulse data, re- 
gardless of amplitude, the wave map evolution has the 
pulse maximum and energy density maximum initially 
moving inward (decreasing r) . For small (subcritical) val- 
ues of A, this inward motion of the maximum proceeds 
for a finite time, after which the maximum "bounces" 
away from the origin and begins to move outward (see 
Fig. |l|). There is a general dispersal of the energy den- 
sity; and for large t, there is very little energy density 
remaining near the origin. 

For large (supercritical) values of A, the behavior of 
the evolving wave map is qualitatively the same initially. 
However, rather than bouncing away from the origin, the 
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maxima for supercritical data continue to approach the 
origin (Fig. |l|), with the concentration of energy around 
the origin appearing to grow without bound. As the en- 
ergy density and the gradient of the function x grow very 
large at the origin, the numerical evolution inevitably be- 
comes unable to resolve the gradient, and the solution 
becomes sufficiently non-smooth to cause the numerical 
evolution to stop. If this accumulation is indeed a sin- 
gularity forming, there is no hope for the numerical evo- 
lution to resolve it, being itself of finite resolution. The 
task then is to examine the behavior of the numerical 
solution up to this point. 

Before doing so, we first discuss a couple of standard 
tests of a numerical solution. We let x(r, t) be some solu- 
tion to the (continuum) partial differential equation (|^) 
and let Xhi'r, t) be the solution to a discrete form of that 
equation, for corresponding initial data, on a grid spacing 
h = Ar. The hope is that, as the grid spacing Ar gets 
smaller, the solutions to the discrete equation generated 
by the evolution code converge to the solutions of the 
PDE, xl^'ji) ~^ x{''':t)- Because in general the explicit 
solutions to the PDE are unknown, we instead consider 
a series of numerical solutions on grids of increasing res- 
olution, say X4h, X2h, Xh- If these are to converge to the 
PDE solution, then they must converge themselves. To 
examine this convergence, we define a convergence factor 
(Q) as follows 



Q= 



IX4h— X2h|2 
\X2h-Xh\2 ' 



(9) 



where the norms are the I2 norm. For these solutions to 
converge, the difference between solutions for increasing 
resolution must decrease and hence Q must be greater 
than one. For second-order schemes, Q is expected to 
be 4. 

Another common test of numerical accuracy focuses 
on the degree to which energy is conserved by the 
numerical evolution. The evolution governed by the 
PDE (^) does conserve energy; the question is whether 
this remains true for the numerical evolution. Letting 
Enum{t) denote the energy calculated from the numer- 
ical solution at time t (on the finite grid), and setting 
g„„„x(t)-g„„„(o) 



A{t) = In 



we monitor A{t) for differ- 



ent choices of grid spacing. The expectation is that A(t) 
should decrease with increasing resolution; if we observe 
this, our confidence in the accuracy of our numerical so- 
lution is enhanced. 
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FIG. 2. Results of a super-critical evolution for an initially 
ingoing Gaussian pulse (^1 = 2, 7?o = 10, S = 2.3, i?max ~ 30). 
The results are shown for increasing resolutions n = 2* (solid. 



(short dash, green), 
(dot-short dash, red). 



black), n = 2-^ (dot, blue), n = 2 
n = 2^^ (long dash, yellow), n = 2^ 

n — 2^^ (dot-long dash, magenta), n = 2^** (short dash-long 
dash, cyan), n = 2^^ (solid, black), and n — 2^® (dot, blue), 
where h = Rmn^/n. The top frame shows the rapid growth 
of x'(0,t) near the time of the blow-up {t ~ 8). The mid- 
dle frame shows the convergence factor (defined in Eq. (^)). 
Factors greater than one indicate convergence. The bottom 
frame shows the change in energy with respect to the initial 
energy. As the resolution increases, so does the level of energy 
conservation. 

In Fig. 1^, we show the evolution in time of three 
quantities-lnx'(0, t), Q(i), and A(i)-for numerical runs 
of supercritical ingoing Gaussian pulse data, done with 
nine different grid spacings. In the top frame, we show 
the behavior of the derivative of x the origin as a func- 
tion of time. The figure shows that as the pulse travels 
inward, the derivative increases. Until just before t — 
all the resolutions show the same behavior as would be 
expected for a convergent evolution. However, near the 
blow-up time, the solutions diverge with higher resolu- 
tions providing a larger derivative. The convergence fac- 
tor Q{t) is shown in the middle frame; it likewise shows 
second-order convergence up to times close to the blow- 
up time. In the bottom frame, the change in energy A{t) 
is shown. We see that as the resolution is increased, en- 
ergy conservation improves. 

What does this tell us about singularity formulation in 
wave maps evolved from (supercritical) ingoing Gaussian 
pulse data? We first argue that these results are con- 
sistent with what would be expected for such formation. 
As the singularity forms, higher and higher frequency 
components become important, and they are represented 
numerically only if one uses higher and higher grid res- 



3 



olutions. Hence, the behavior of the derivative of x ^-s 
the resolution improves would be expected to show larger 
and larger gradients, as seen in Figure 2a. Next, we note 
that the formation of a singularity should not hinder con- 
vergence except quite near the formation time, as is seen 
in Figure 2b.. Finally, energy conservation should be fine 
until the high frequency components play their role, as 
we see in Figure 2c. Hence, the results observed appear 
to be consistent with a singularity forming near t = 8. 

This does not guarantee that a singularity forms in 
these wave maps. There are other effects that might pro- 
duce the apparently unbounded growth of the derivative 
of X of the energy density near the origin in these 
numerical simulations. For example, perhaps some un- 
physical, unstable mode grows because of the particulars 
of our chosen evolution scheme. We believe that this is 
not the case, for a number of reasons. First, the presence 
of such a mode would likely cause much larger growth 
in A(i) than we see. Second, such modes would have to 
be excited only after some time (roughly independent of 
resolution) and only for families of sufRciently large en- 
ergy. This is not consistent with our observations. Third, 
the excitation of this sort of instability would almost cer- 
tainly depend critically on the precise finite difference 
scheme. Because Bizoh and his collaborators ob- 
serve similar behavior, using a different numerical evolu- 
tion scheme, this does not appear to be the case. Thus we 
believe it very unlikely that the effects we are seeing are 
the result of a numerically unstable nonsingular mode. 

Another situation in which one might numerically ob- 
serve the formation of singularities that do not in fact 
evolve analytically from the corresponding data is if the 
continuum PDE solution is regularized by high frequency 
components which cannot be seen by the finite grid res- 
olutions we use. The rather strong convergence behav- 
ior we see in our numerical solutions leads us to believe 
that this is not happening. We note in particular that 
such unresolved components would have to be separated 
in frequency space from the nontrivial low frequency 
components by a substantial margin, with a large range 
of dynamically irrelevant frequencies separating the two 
regimes. This seems to be very unlikely. 
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FIG. 3. Near critical evolutions approach the static solu- 
tion, Eq. (0). Nine frames equally spaced in time are shown 
for both suh- and SMper-critical evolutions. The solutions are 
indistinguishable in the graph at these times and are shown 
with dots. After t — 14.24, the two solutions have quite differ- 
ent fates, but both approach the form of the static solution. In 
the final frame, the static solutions x(r) ~ — 2 arctan (lllGr) 
is shown. 



IV. CONCLUSION 

The numerical studies we present here very strongly 
support the contention, previously conjectured by Shatah 
and Struwe ||^, that smooth initial data for wave maps 
from 2-1-1 Minkowski spacetime into the round two sphere 
can develop singularities (with unbounded derivatives) in 
finite time. As we note, there are many ways in which 
the numerical exploration of possible singularity forma- 
tion might produce misleading indications. However, we 
believe that as a consequence of the numerical tests we 
have carried out, together with those done independently 
by Bizoh and his collaborators , the formation of sin- 
gularities is the most likely conclusion. 

There is much more one would like to know about these 
spatially equivariant wave maps, as well as about those 
without such symmetry. One would like to know, for 
example, if the solutions assume any universal form as 
one approaches the singularity. Our work (see Fig. |^) 
supports the results of Bizoh et al which indicate 
that indeed the family of static spherically equivariant 
wave maps 

X(r) = ±2arctan(Ar) (10) 

does serve as a sort of universal model for singularity 
formation. This needs to be studied further. One would 
also very much like to understand the behavior of the 
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wave maps which evolve from initial data near the tran- 
sition from subcritical to supercritical data. The recent 
numerical work of Bizon et al [|lO| suggests that the 
static solutions (^0|) play a central role in the evolution 
of the transitional wave maps as well in that of supercrit- 
ical ones; however, this issue needs further investigation. 
(Note the absence of any self-similar solutions to the 24-1 
wave map equations; for 3-1-1 wave maps, such solutions 
play a key role in behavior of solutions evolving from 
critical or near critical data). 
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